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Abstract. We present an analytic approach to study concurrent influence of quenched non-magnetic site- 
dilution and finiteness of the lattice on the 2D XY model. Two significant deeply connected features of this 
spin model are: a special type of ordering (quasi-long-range order) below a certain temperature and a size- 
dependent mean value of magnetisation in the low-temperature phase that goes to zero (according to the 
Mermin- Wagner-Hohenberg theorem) in the thermodynamic limit. We focus our attention on the asymp- 
totic behaviour of the spin-spin correlation function and the probability distribution of magnetisation. The 
analytic approach is based on the spin-wave approximation valid for the low-temperature regime and an 
expansion in the parameters which characterise the deviation from completely homogeneous configuration 
of impurities. We further support the analytic considerations by Monte Carlo simulations performed for 
different concentrations of impurities and compare analytic and MC results. We present as the main quan- 
titative result of the work the exponent of the spin-spin correlation function power law decay. It is non 
universal depending not only on temperature as in the pure model but also on concentration of magnetic 
sites. This exponent characterises also the vanishing of magnetisation with increasing lattice size. 
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PACS. 05.50.+q Lattice theory and statistics (Ising, Potts, etc.) - 64.60.Fr Equilibrium properties near 
critical points, critical exponents - 75.10.Hk Classical spin models 

1 Introduction The spin- wave approximation (SWA) applied by F. Weg- 

ncr [4] to analyse the 2D XY model leads to a result very 

The quasi-long-range ordering (QLRO) is a special feature c i ose to recent Monte Carlo computations in the region of 

of a number of important many-particle systems including \ ow enough temperatures. In particular, the presence of a 

two-dimensional solids, magnets, Bose fluids, liquid crys- special type of ordering - the QLRO - manifests itself in 

tals PEE] . As it is known by now the 2D XY model serves the power law decay of the spin-spin correlation function: 
as an archetype capturing special features of QLRO in 

these systems. Here we will focus on this particular model (S r S r +R,) ~ i?~''° B (2) 
keeping in mind that the results can be generalised for 

some other similar models. The regular model, described wher <3 R is the distance between the spins. The exponent 

by the Hamiltonian 7 f cg given by the SWA is non-universal: 

j 7f cg = kT/(2nJ) . (3) 

-^rcg 2 r )(^r^r'+^r^r') ' (-0 ^he detailed description of properties of the model 

given by V.L. Berezinskii [5], and J. M. Kosterlitz and 

has been investigated in great detail, and although most D. J. Thouless (BKT) [6|7j is based on the hypothesis 

of its properties are known, no exact solution was found, that certain local spin configurations, named topological 

In ([1]) r and r' span sites of a two-dimensional lattice, defects, are responsible for the QLRO and the behaviour 

J(r) is the nearest neighbours interaction potential, , of the system near the transition to the QLRO phase (the 

S% are the components of a classical "spin" S r , the coef- BKT transition) at the temperature Tbkt- The intuitive 

ficicnt 1/2 stands to prevent double count of each bond, analogy with the transition in electrolytes was used in 
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those works. A further analytical basement for this ap- 
proach can be found in the work of J. Villain [8j. 

There are two important aspects which differ the ideal 
2D XY model from the systems that can be met in na- 
ture: real physical systems are always finite and possess 
structural defects. The finite size effect already has been 
widely explored in the pure (undiluted) XY model. The 
interest of this question is that in any 2D XY system of 
finite size, magnetisation is non-vanishing [9] and goes to 
zero only when the lattice becomes infinite as it should 
be according to the Mcrmin-Wagncr-Hohcnbcrg theorem 

mm- 

The decay of the magnetisation has a power law form 
(below the BKT transition temperature) that can be eas- 
ily found in the SWA PTT5] : 

(m) - N- 1 ^ , (4) 

with the same rj les , Eq.([3|), that stands in the correlation 
function @. Here, N is the total number of spins. The 
same result can be obtained also from the finite size scaling 
(see e.g. Ref. Q3] in a similar context.). 

The recent works of S. T. Bramwell et al. [13115116) 
give deep analysis of the magnetisation probability distri- 
bution which they claim to be non-Gaussian and of uni- 
versal form, independent of both system size and critical 
exponent ?7 rcg . 

Structural disorder as site- or bond-dilution deserves 
much attention since it moves an ideal model closer to- 
wards true physical systems which can be found in nature. 
However the number of works dedicated to this aspect in 
the 2D XY model is not mirrored in the great impor- 
tance of the topic. Harris criterion [17] implies that energy- 
coupled disorder has no effect on the universal properties 
(e.g. the critical exponents) at the transition temperature. 
The BKT universality class (and in particular the cele- 
brated t](Tbkt) = j) is thus unchanged by the introduc- 
tion of quenched disorder, but one can expect highly non- 
trivial dependence of the low-temperature characteristics, 
like the spin-spin correlation function, on the concentra- 
tion of spin- vacancies [21j . 

A non-magnetic site can change the interaction be- 
tween topological defects which are responsible as it was 
mentioned above for the QLRO |18|19|20] , The charac- 
ter of this influence is not completely clear up to now, for 
example, the question: when the QLRO disappears as the 
concentration of vacant sites increases, has got different 
answers [20121122] , As it appears now, the most convincing 
scenario is that the QLRO remains up to concentrations 
very close to the percolation threshold |21I22| . However we 
do not touch this question focusing mostly on the region 
far from the percolation threshold. 

In this paper we investigate the concurrent influence 
of quenched site-dilution and finite size of the lattice on 
the properties of the 2D XY model. These two modifi- 
cations together present a nice approach to investigation 
of real physical systems. To quantify the disorder-induced 
changes in the QLRO phase we pay attention to the spin- 
spin correlation function exponent of the 2D XY model 
with quenched site-dilution, r/ dl1 . It describes not only the 



decay of the correlation function with the distance but 
also the vanishing of the magnetisation in a finite system 
with increase of the lattice size and the divergence of the 
susceptibility in the same limit. The analytic approach we 
use here relies on the SWA and a perturbation expansion 
and is verified by MC simulations. 

Also we will perform Monte Carlo simulations for a 
wide range of 2D A"F-spin systems of different sizes at 
different temperatures and with different concentrations 
of impurities. Systems explored in computer experiments 
are always finite, thus they possess a non-vanishing mean 
value of magnetisation. The instantaneous magnetisation, 
which is the scalar value of the total sum of the spins di- 
vided by the number of sites, measured in a given state 
from the thermodynamical ensemble of states of the sys- 
tem is distributed with a certain law. The form of this 
distribution is the point of our interest as well. 

The structure of the paper is the following: In the sec- 
ond section we give a description of the model and calcu- 
late analytically the spin-spin correlation function combin- 
ing the SWA and perturbation expansion, we support the 
result by MC simulations. In Section 3 more details about 
MC simulations can be found. The results are presented 
in the form of ring functions and probability distribution 
functions of magnetisation. We analyse the plots and add 
an analytic calculation of the moments of magnetisation. 
We discuss the analytic and MC results and sketch the 
plans of future work in Conclusions and give two appen- 
dices with technical details of the calculations. 



2 The spin-spin correlation function 

In this section we give description of the diluted 2D XY 
model and explain the expansion applied to analyse the 
asymptotic behaviour of the spin-spin correlation func- 
tion in the low temperature limit. The comparison with 
our Monte Carlo results is added to support the analytic 
approach. 

2.1 The model 

The regular 2D XY model ([5]) is equally described in the 
angle variables £? r 's that are the angles between the spins 
and a certain fixed direction by the Hamiltonian 

ffro g -^E - r ') COS (^ -8*), (5) 

r r' 

since S^S*, + S^S^, = cos(# r — 9 r >) for unit length spins. 
All the notations are the same as in ([I]). 

We define a set of occupation numbers c r 's that intro- 
duce disorder into the lattice: 

{1, if the site r has a spin; 
0, if the site r is empty. ^ ' 

The Hamiltonian modified with these numbers, 
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where 



H = -- ^ E ~ F ') C ° S ^ r ~ MCrCr 



(7) 



will describe the model on a lattice with dilution. Setting 
a certain sequence of numbers {c r } we are able to real- 
ize any configuration of lattice dilution with any desirable 
concentration of magnetic sites c = c7- We are interested 
in some thermodynamical quantities which are dependent 
in this case on the configuration of impurities. To obtain 
observable values we will average the quantities of interest 
over all the possible configurations of non-magnetic sites; 
this is referred to as quenched disorder in the literature 
as in contrast to annealed disorder when magnetic and 
non-magnetic sites are in equilibrium and the configura- 
tional averaging has to be made already in the partition 
function [22] . We denote the configurational averaging as 



(...) 



II E [cSi- Ct ,o + (l-c)5 Crt0 ](...) (8) 



r Cr=0,l 



(Si j are Kroneker deltas), here and below index r in sums 
and products spans all sites of a 2D square lattice. 

Since we restrict ourselves to the low temperature phase 
of the model we assume that the directions of spins on 
neighboring sites do not differ essentially This allows us 
to pass to the SWA replacing cos (9 r — 6 T i ) in the Hamilto- 
nian with a quadratic form 1 — | (6 T — 6 T i) . All the main 
features of the model are preserved in the low temperature 
limit in the Hamiltonian 



(9) 



-«qr 



(c r - c) , 



5k.k' 

7k = 2 



7k 

— cos k x a — cos k y a 



(14) 

(15) 
(16) 

The last relation is true for a square lattice with spacing a 
and the nearest neighbours interaction of strength J. It is 
important to stress that the first term in the Hamiltonian 
can be regarded as the SWA Hamiltonian of the model on 
a pure (undiluted) lattice with a rcnormaliscd coupling. 
We can use it and write thermodynamical averaging with 
respect to the Gibbs distribution as 



e -(3(H p + H p2 )\ /-f)(H p +H p2 \ 



(17) 



where the notation (...)* is used for thermodynamical av- 
eraging with the Hamiltonian of the pure system: 



Tr 



-c 2 /3H r{ 



Tr e" c2 ^< 



(18) 



Tr (...) 




+00 



d9i 



-foe 



dot (...) , 



where Q% = $19 T , 9^ = 3# r , and in order to keep the same 
number N of variables the product has to be taken over 
a half of the 1st Brillouin zone which we have denoted as 
B/2 . It was possible to extend the integration region to 
(—oo, +00) because of the Gaussian form of the Boltzman 
factor that stands in the integrals. 



where the first term in the expression can be regarded just 
as a shift in the energy scale, from now on we denote the 
second term by H for simplicity. 

Using Fourier transformation of the variables: 



2.2 The expansion in {p q } 



^E* 



e'k' 9\r. 9\r 



J ( r ) = jf E e *Mq)> Kq) = E e "' qr J ( r ) 



Let us note, that the transformation (fTTjl of the Hamil- 
tonian (J9j) is exact, although it looks like a perturbation 
expansion in p. In the forthcoming calculations in order 
\ g-ikr^ to perform configurational averaging we expand any ther- 

^ ' modynamical quantity of interest (F({pq\)) in terms of 

functional variables pk's, Eq.(frj~ 



where k runs over the 1st Brillouin zone, one arrives at: 



H = c 2 H [cg + H p + H p 2 , 



with 



7kSk,k'Pk+k'tfkPk' 



(11) 



(12) 



(%})) = (f({o})> + EM k K (19) 

k 

+E/2( k > k 'Wk'+ y, / 3 (k,k',k'> k 



k.k' 



k.k'.k" 



k,k' 



k,k',q 

X [P-k-k'-qPq - P-k-qP-k'+q] #k#k' (13) 



Since p^s characterize the deviation from the com- 
pletely homogeneous disorder in the Hamiltonian they can 
be considered as parameters of perturbation. Note that a 
power of p corresponds to the number of sums over k in 
([i~9|) . A classification of the perturbation theory series with 
respect to the number of sums over k corresponds to the 
expansion in the ratio of the volume of effective interac- 
tion to the elementary cell volume [2"4] . Taking this ratio 
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to be small means that it is valid for the short-range inter- 
acting systems, which holds for our problem. As far as we 
don't make any assumption about weakness of disorder, 
we may expect that accordance of results of this expansion 
with the MC simulations should not be very sensitive to 
the value of dilution (1 — c) (but of course still far enough 
from the percolation threshold where the whole approach 
fails). 

In the calculations presented below we limit ourselves 
to the third order term in the expansion. Then it is not dif- 
ficult to perform averaging over configurations of disorder 
using the equalities: 



Pq. = , 

PqPqT = c(l-c) — <J q+q ,, , (20) 

PqPq'Pq" = c( 1 - 3c + 2<?) <S q+q < +q » ,0 

which can be obtained easily from (JSj> . 



2.3 The asymptotic behaviour of the spin-spin 
correlation function 

The spin-spin correlation function of the diluted 2D XY 
model, 



G 2 (R) = (c r c r+R cos(0 r+R -Or)) , (21) 
can be written in the Fourier variables (1101) as 



with 



c r c r+R cos -L [VIK + r£6t) 



?/£ = coskr — cosk(r + R) , 
77^ = — (sinkr — sink(r + R)) 



(22) 



Writing the expansion (|19[) and applying the equalities 
(|20p we arrive (see appendix A) at the next expression: 



G 2 (R) = c 2 (cos(0 r+R -# r )), 

1 -c 1 1 ^ 



1 - — JjW^ 9 ^' 9 * M ~ 



1 - 3c + 2c 2 l/2^ 
~ A ^ V* k 



N 3 ^ 

k,k',k 



\ ^ sin 
<?-k,k'S'k',k"fl , k",k — ~ 



(23) 



Since Eq. (f23"l) is already configurationally averaged it 
docs not contain the p's anymore, so the correspondence 
with the orders of the expansion Q19p is not obvious. Let 
us explain the origin of each term. The unity corresponds 



to the zeroth-ordcr in p, the first-order term is identically 
vanishing as follows from (|20[) . and the second- and third- 
order terms in p can be distinguished by their coefficients 
that are clear from (f20|) . The second-order term contains 
two sums according to the expansion (jT^J) ■ At the same 
time the third-order term contains, except the triple sum, 
a sum over one k: the summation over two remaining k's 
was possible to carry out explicitly in this particular case. 

For our purpose it is enough to get the leading asymp- 
totics of the sums that stand in the expression ([23"]) when 
N — > 00 and R — > 00 (see Appendix B): 



N ^ 



sin 2 _ 1 , fj 



k 

1 

iV 2 



y^gk.k'ffk'.k 



const + — In — , 
2tt a ' 



~s const' 



7k 



k.k' 



73 .9-k.k'3k',k"gk",k 



N 3 



2tt 
const' 



ln^ , 

n 1 



im 2 ^ _ i „ 0.27, R 



k.k'.k" 



2- 



ln- 



Inserting these expressions in (|2"3")> it is possible to write 
the pair correlation function for small enough tempera- 
tures in the power law form: 



G 2 {R) « c 2 (R/aY 



(24) 



Reminding the spin-spin correlation function exponent of 
the pure system, rj les , given in the SWA by Eq.Q, we 
write 



n 



an 



1 



0.73- 



1 



2.27 



(l-3c+2c 2 



(25) 

In fact, as it can be seen from Appendix B, this result 
is true not only in the thermodynamic limit but also for 
a system of finite large enough size N. The first term in 
the brackets, 1/c 2 , corresponds to the zeroth order of the 
p-expansion, the first-order term is identically vanishing 
as was already noted before, the second and third terms 
in the brackets correspond to the second- and third-order 
terms of the p-expansion respectively. In the next subsec- 
tion we evaluate formula (|25[) and compare this result with 
the MC experiments. 



2.4 Comparison with the Monte Carlo results for the 
exponent of the spin-spin correlation function 

In order to check Eq. (|25j) . we have performed simulations 
of 2D XY-spins using Wolff cluster Monte Carlo algo- 
rithm [25] . We only mention here the main features of the 
simulations used in order to obtain the exponent ?y dl1 . We 
discard typically 10 5 sweeps for thermalization, and the 
measurements are performed with typically 10 5 produc- 
tion sweeps. Averages over disorder are performed using 
typically 10 3 samples. There is no need of a better statis- 
tics, since we are far from the BKT point (in the vicinity 
of the deconfining transition, the presence of many topo- 
logical defects is an obstacle to thermalization as it can be 



O. Kapikranian, B. Berche and Yu. Holovatch: The 2D XY model on a finite lattice with structural disorder 



5 



shown empirically by the analysis of autocorrelation time 
(see e.g. Ref. [26])). The boundary conditions are chosen 
periodic and the critical exponent rj(T) of the correlation 
function is measured indirectly through the finite-size scal- 
ing behaviour of the magnetisation 



M T (L) ~ L~ x '^ , 



MT) = \v{T), (26) 



where the last scaling relation holds in two dimensions 
(L = \/~N is the linear size of the lattice). 

In Fig.[T] we compare the ratio r] AA /rf eg evaluated an- 
alytically by keeping terms from the 1st to 3rd order in 
p-expansions with the MC data. One can see that up to 
the third order the analytic curves approach step by step 
the MC data: the Oth order seems to be a rough approx- 
imation, the 2nd order curve lies closer to the MC data, 
but still much below, and the third order curve seems to 
fit the MC results better, although it doesn't give perfect 
accordance. Of course this fact does not allow to conclude 
in favour of a similar agreement for higher orders, and it 
is possible that the expansion will not show any conver- 
gence at all. On the other hand, having more perturbation 
theory contributions at hand one can attempt to apply a 
resummation technique to improve its convergence, sim- 
ilarly as it is commonly done analyzing field-theoretical 
expansions [57J . 






Oth order 




2nd order 




3rd order 


o 


MC 0.04 


□ 


MC 0.08 


o 


MC0.12 



Fig. 1. The comparison between the MC data and the analytic 
results for the ratio ?? dll /?7 reg as a function of concentration of 
occupied sites c obtained in different orders of the p-expansion. 



Anyway, comparing outcomes of our analytical and 
MC treatments presented in FigQ]one arrives at the con- 
clusion about good agreement within concentrations from 
c = 0.75 to c = 1 at least up to the third order of the 
perturbation expansion. 

Let us proceed further investigating magnetization and 
its distribution in a finite-size system. 



3 The magnetisation probability distribution 

In this section we obtain and discuss the probability distri- 
bution function of magnetisation obtained in Monte Carlo 



simulations of a two-dimensional XY-spin system, per- 
formed for different sizes of the lattice, temperatures and 
concentrations of impurities. We support the MC analysis 
with an analytic treatment of the magnetisation probabil- 
ity distribution in a model of finite size using the same 
approach as for the spin-spin correlation function. 



3.1 The probability distribution functions 



In subsection 12.41 we obtained the spin-spin correlation 
function decay exponent by use of the finite size scaling re- 
lation ([26]) for the magnetisation measured in MC simula- 
tions. At the same time the probability distribution func- 
tion (PDF) of magnetisation itself deserves much atten- 
tion since it appears to be of non-trivial form. As we men- 
tioned in Introduction it is known that it is non-Gaussian 
in the pure 2D XY model. In the case of structural dis- 
order we can expect also dependence on concentration of 
dilution. 

We define instantaneous magnetisation as the scalar 
value of the total sum of the spins divided by the number 
of sites, 



1 

N 



(27) 



measured in a given state from the thermodynamical en- 
semble of states of the system with a fixed configuration 
of structural disorder. The probability to find the system 
in a state with magnetisation to, P con f (to), considered as 
a function of to is called the probability distribution func- 
tion (PDF) of magnetisation or just the magnetisation 
probability distribution. 

The thermodynamical mean value of magnetisation de- 
fined through the usual procedure of thermodynamical av- 
eraging with the Hamiltonian of the system H: 



Tr (j 



' Tre-i 3H ' 
can be written then in terms of the PDF as 



(m) = / mP con{ (m)dm . 
Jo 

We define also the pth moment of magnetisation as 



mP P coa f(m)dm 



(28) 



(29) 



(30) 



In this place it is important to emphasize that the mag- 
netization to, which stands in the integral in (|29[) . differs 
from that defined by Eq. (|2~7} . since in ([29]) to just plays the 
role of the integration variable which doesn't depend on 
the microscopic state of the system. The PDF, P con f (w), 
is already a thermodynamic characteristic depending only 
on a macroscopic state of the system. It can be seen from 
the well known property of probability distribution func- 
tions that a PDF is defined uniquely by its moments [29] : 
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r°° dx ^> (-ix) p 
P C oAm)= ^J2^L Mp . (31) 

The moments are thcrmodynamically averaged quantities, 
so it follows that the PDF of magnetization is a thermo- 
dynamic quantity too and depends on m only as on a 
parameter, however to find an analytic expression for the 
PDF is not a trivial task even for the pure model. 

Thus the mean magnetisations defined by Eq. (|28|) and 
Eq. (j2"5)) being the same are obtained by different proce- 
dures. 

Since we investigate observable quantities here we should 
look at the configurationally averaged values of magneti- 
sation and its moments. In terms of the PDF it means 
that P C onf (tO must be averaged over the whole range of 
possible configurations of impurities with a fixed concen- 
tration. Then the mean magnetisation and its moments 
can be written as: 

(m) = / raP{ra)dm (32) 
Jo 

and 

TTp = ( m p P{m)dm , (33) 
Jo 

where P(m) = P C onf (jn) is the configurationally averaged 
PDF. 

The PDF of magnetization is very suitable for further 
analysis of results obtained in Monte Carlo simulations. 
In Fig. [2] we illustrate the procedure of configurational 
averaging of MC data. From different curves of P C onf {rn) 
obtained for different realizations of dilution we draw one 
averaged curve which is P(m). 




0.9 0.95 

m 



Fig. 2. The probability distributions of magnetisation for 
twenty different realizations for a system of size L = 16 at 
dilution c = 0.95 at a temperature ksT/J = 0.1. The thick 
line is the average probability distribution, still very bumpy 
with so few configurations. 



°- 07 | 




0.2 0.4 0.6 0.8 1 



m 

Fig. 3. The average distributions over 10 3 samples at temper- 
atures ksT/J = 0.1, 0.5 and 0.9 for systems of increasing sizes 
L = 16, 32, 64. 

In the pure 2D XY model two of the main features of 
the PDF are that the form of the distribution at fixed tem- 
perature is universal, i. c. it does not depend on the size 
of the system, and it is non-Gaussian. These two state- 
ments have been derived analytically and verified by MC 
simulations [13I16|15] , 

Fig. [3] illustrates the MC results for size dependence 
of the PDF of magnetisation in a diluted 2D XY-spin 
system with concentration of spins c = 0.95 at three dif- 
ferent temperatures. We see that at fixed temperature the 
mean magnetisation becomes smaller for bigger lattices as 
it should be according to Eq. ((4]) in the case of a pure sys- 
tem. From the curves for low temperatures it is clear that 
the form of the distributions is non-Gaussian like in the 
pure model. What is more interesting is that the form of 
the distributions is noticeably different for different sizes. 
It seems to be in contradiction with results for the pure 
model [13I16115| . 

A suitable parameter that can characterize the form 

of a PDF is the variance: a = \J M2 — M\ . It has been 
prooved that in the pure 2D XY model it is independent 
of system size. Here we see in Fig|3]a different qualitative 
behaviour, the variance, which is proportional to width 
and flatness of the distribution, grows as the size of the 
lattice decreases. 

Since we simulated a diluted system, this must be the 
result of non-magnetic impurities influence. It calls for an 
analytic explanation of this dependence in Subsection 3.3. 



3.2 The ring functions 

Another way to display the magnetisation probability dis- 
tribution observed in MC simulations is to draw a ring 
function which is denned in the next way. A ring function 
is obtained when one plots the successive values of the 
magnetisation (for each Monte Carlo step) in the plane 
(m x ,m y ) where m x and m y are the two components of 
the magnetisation (see Fig|J). In fact it contains the same 
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information as a PDF. the difference is that in a ring func- 
tion we can see the distribution of the vector of magneti- 
sation, its form however is symmetric around the point 
rn = 0, i. e. there is no preferable direction (this is the sig- 
nature of rotational invariancc of the QLRO phase), and 
the probability, P(m x ,m y ), is mirrored in density of the 
points. 

Since the algorithm used here is a cluster algorithm 
specially dedicated to this type of spin systems, the suc- 
cessive spots are essentially uncorrclated. This is a very 
different situation with a Metropolis algorithm where the 
successive spots would be correlated (see Ref. [13j). 

We are interested in the temperature dependence of a 
ring function of magnetisation in a diluted 2D XY-spin 
system with fixed size and concentration of impurities. 
For this purpose ring functions for a system of size L = 
16 at dilution c = 0.95 for three different temperatures, 
ksT/J = 0.1, 0.5 and 0.9, are shown in Fig.|U The outer 
ring functions (color on-line) represent the pure system 
of the same size and at the same temperatures. The first 
feature that one can notice is that the radius of the rings 
of higher density of the "diluted" ring functions is smaller 
for all temperatures than that of the corresponding "pure" 
ring functions. This is due to the fact that we consider the 
magnetisation per site taking all sites into account and in 
a system with impurities there are missing spins. 

In which concerns the temperature dependence of the 
ring function, it is visible that the mean magnetisation 
tends to zero as the temperature increases in both pure 
and diluted system. When for the pure 2D XY model 
this behaviour follows from Eq.((4]) (since ?7 rcg ~ fcsT), 
it must be verified analytically for the case of a model 
with site-dilution, what is done in Subsection 3.3. An- 
other feature of the temperature-behaviour which can be 
noticed in FigfJ]is that the high-density region is wider for 
higher temperatures and approaches to a delta-function as 
the temperature goes to zero, this feature is well known 
from MC and analytic investigations of the pure 2D XY 
model piT5j . 

Except the difference in position of the density peaks 
of the ring functions caused by the difference in the total 
numbers of spins, there are no qualitative differences in the 
form of the "diluted" and "pure" ring functions obtained 
at the same temperature for the value of concentration of 
magnetic sites c = 0.95. Their widths are approximately 
the same and both distributions show clear non-Gaussian 
character for higher temperatures which lies in the visible 
fact that more points arc situated in the inner region of the 
ring functions than out of the rings of the highest density. 

The discussion above concerns relatively weak dilution 
with c = 0.95. The situation appears quite different when 
the dilution becomes stronger. In Fig. we show a ring 
function for a system of the same size as in the previ- 
ous figure at temperature kgT/J = 0.1 but with much 
stronger dilution, c = 0.70. One can see now that in a sys- 
tem with sufficient number of non-magnetic impurities the 
high-density region of the ring function is much wider than 
in a pure system at the same temperature. We can con- 
clude thus that the variance that characterizes the width 
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Fig. 4. Ring functions for a system of size L = 16 at dilution 
c = 0.95 for temperatures (from top to bottom) ksT/J = 0.1, 
0.5, and 0.9. The outer ring function (color on-line) represents 
the pure system. 



of the distribution must be dependent on concentration of 
dilution in a 2D XY model with structural disorder. This 
non-trivial observation calls for an analytic support which 
is the point of the next subsection. 
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Fig. 5. A ring function for a system of size L — 16 at dilution 
c = 0.70 for temperature fcsT/ J = 0.1. The outer ring function 
represents the pure system. 



3.3 Analytic calculation of magnetisation and its 
moments 

In this subsection wc present an analytic analysis of the 
PDF of magnetisation in the 2D XY model in the case 
of structural disorder. We use here the same scheme that 
was presented in Section [5] to obtain the moments of mag- 
netisation which define the magnetisation probability dis- 
tribution through Eq. (f3Tj) . 

It is convenient for analytic calculation to rewrite the 
instantaneous magnetization (|27p as 



cos 



(34) 



where 8 = jj^2 r & r is the algebraic mean of the angle 
variables. 

Since we are interested in observable quantities wc 
must consider configurationally averaged moments of mag- 
netization: 



Mr, 



(™ p ) = ( Cri ■ ■ ■ Cr p C0S7/;r 



cos 



where we denoted ipr = T — 9. Note here that the 1st 
moment of magnetisation, 



1 = — ^2 ( Cr C0S7 M = (co cos ipo) , 



M 



is just the mean magnetisation (m) by definition. 

Passing from the product of cosines to a sum we write 
the (n + l)th moment as 



+ i 



(35) 



The expression that stands in the sums in (|35[) can be 
written in Fourier variables as: 



(c c ri • • • c r „ cos(?Ao + aiipn ha„^ r J) 



c c ri • • • c r „ cos ^j- £ k (f£0£ + j£0£) 



with 



rjy = 1 + a.\ cos kri + ■ ■ ■ + a n cos kr n , 

= — («i sinkri + • • ■ + a n sinkr n ) . (36) 

Substituting in ([35)) the result for the above expression 
obtained in Appendix A in the third order approximation 
in the p-expansion, Eq. lfSD)) . we get 



„n+i 



M n+1 



1 - 3c + 2c 2 



2 n N 

ri , . ..r n a.i— ±1 

1 (l-c 1 



7 



9k', k 



' \ k.k' 



27V Z 47V3 Z 

k k,k',k' 



5k,k'5k',k"5k",k 



X f n + 1 + 2 J2i<j a i a 3 cos k ( r i _ r j 



7k 



(37) 



with .gk, k' given by Eq. (fT5)) . Then, using the result for 
a quantity of the type ^cos ^= ^ k fak^k + V^L)) ^ from 
Appendix A, Eq. (|45[) . we write in the low-temperature 



limit 

(COS^O + aiVW ^nVO), 



(38) 



ac4jn J: ( n + 1 + 2 Si<j cosk ( r i - r j) 

k^0 



Substituting this expression in (|3"Tj) we are able to sum 
over all a^s by help of the obvious equalities: 



^2 on = , of = 1 . 



a i= ±l 

Then the (n + l)th moment of magnetisation reads as 
n + 1 / 1 1 



Mn+l 
n+1 = C 

1-C 1 



1 - 



/3 J \ 4c 2 N ^ ^ k 

M \ k^O 



(39) 



4c 3 N 2 



X].9k,k'.9k',k^ 



k,k' 



1 - 3c + 2c 2 
2c 1 



TV 51 7k 2A 3 Z 

v k k,k',k" 



fl , -k,k'ffk',k"3k",k^ 



In the limit N — > oo we find for the sums in ([3T))) (see 
Appendix B): 



2 ™^r I] Z ^ c ° c ^ • • • cos ^° + «i^i)> ■ ^ z 4 ~ const + 27 lniV ' 

ri ,. . .r n ai=±l 
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AT2 E E ■' 



E 5 ~ k ' k ' 



i , 0.73, Ar 

N 2 ^ ^ ffk,k'5k',k^ « const + — In JV, 

k k' 

i ,, 0.27 

.9k'.k"3k".k— ~ const - 



N 3 



In TV . 



k,k',k" 



For small enough temperatures it is now possible to write 
the pth moment of magnetisation in the form: 



M„ ps c^AT S" d " 



(40) 



with the exponent ?y dl1 given by Eg. (|2"5|). The equality (gDJ) 
can be rewritten as 



M p ps Mi 



(41) 



As far as all higher moments M p can be trivially expressed 
in terms of Mi this relation implies absence of multifrac- 
tality and it differs from that for the pure model in [TS] : 




1 



n(n — 1) ^ 1 
(/3J) 2 167V 2 ^ ~ 



7 2 



since we neglected all the terms in the expansion contain- 
ing powers of 1/((3J) higher than one. In fact the relation 
([4Tj) corresponds to a delta distribution of the probability 
of magnetisation (when the variance is equal to zero) that 
is close to the truth only for very low temperatures. 

At the same time the mean value of magnetisation 
which can be obtained from Eq. (|4"01) in the particular case 
p=h 

Jmj ps ciV~3'? dil , (42) 

recovers the finite-size scaling relation (|2"fJ)) and accords 
well with our MC simulations (see Fig[T]). 



4 Conclusions 



the exponent converges to the MC data with every next 
order. 

We also took into account the finiteness of the lat- 
tice which appears to be negligible for the exponent 7y dl1 
but brings on stage another important property of the 
2D XY model: non-vanishing magnetisation that tends to 
zero with a power law as the size of the lattice increases. 
Monte Carlo simulations of diluted 2D AY-spin systems 
with different sizes, concentrations of dilution and tem- 
peratures, presented in terms of magnetisation probabil- 
ity distribution, show some interesting features that differ 
essentially from the case of the pure model in the same 
conditions. We observed that the variance of the proba- 
bility distribution function of magnetisation, which serves 
as a characteristic of the distribution form, i. e. depends 
on its width and flatness, is a function of temperature, 
system size and concentration of dilution in contrast to 
the pure 2D XY model where it depends on temperature 
only. 

We applied our analytic approach to compute the mo- 
ments of magnetisation that define the probability distri- 
bution function (and its variance) but failed to explain 
the features present in our MC simulations because of the 
roughness of our approximations. At the same time the 
analytic calculations give a good result for the magneti- 
sation itself in the diluted model, we found the power law 
decay with system size, Eq. (|42|) . that accords with the 
finite-size scaling (f2"o]l . 

The convergence of the perturbation expansion applied 
in the analytic treatment is not undoubted but from our 
results we can conclude that up to the third order it gives 
nice accordance with the MC data. Further analysis of 
this question will be a subject of a separate study. An- 
other direction of future work would be to implement the 
same type of perturbation expansion within the Villain 
model [5] and to explore the deconfining transition of the 
diluted model. Here the interesting part of the question 
does not come across the value of the exponent rj, but 
rather in the way the vacancies couple to the unbinding 
mechanism. 



Two important modifications, quenched site-dilution and 
finiteness of the lattice, were brought to the usual 2D XY 
model in order to investigate quasi-long-range ordering, 
which appears in this model at low temperatures, in con- 
ditions closer to that in real many-particle systems present 
in nature. 

We proposed a method of an analytic treatment of 
the 2D XY model with structural disorder based on the 
SWA and a sort of perturbation expansion in the param- 
eter characterising deviation from the pure system with 
the renormalised coupling strength. Computing the per- 
turbation expansion up to the third order we arrived at 
the result for the spin-spin correlation function decay ex- 
ponent ?7 dl1 , Eq. (j25|) . which appears to be non-universal 
and depends besides temperature also on concentration of 
non-magnetic impurities in the system. Our analytic re- 
sult shows nice accordance with MC simulations for a wide 
range of dilution concentrations (Fig[lJ) . We see that up 
to the third order in the expansion the analytic results for 
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Appendix A 

In this appendix we give a detailed calculation of a quan- 
tity of the type: 



(43) 



used while computing the spin-spin correlation function 
and the mean magnetisation of a finite-size system and its 
higher moments. 
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Expanding this quantity in the parameters /9 q 's up to 
the third order we have 



c ri • • • Cr, cos -f- Ek (Tjgfljj. + t£0£ 

= c' (cos^Ek(« + «) 



{Hp cos) t 
(cos) ~ 

(H „ 2 COS) 



or* 



2 / (Hi cos) 



(cos), * ( H P 2 )* ) + 2 ( <cos) t * C^p)* 

2 / (g p fl p 2 CQ S ) t 



(cos) ~ 



(cos) ^ 



(H P H P .) 
-{H p ).(Hf) 



/? 3 // (gp)»(gp 2 cos) t 



(cos) ^ 
z( (H P )l(H p cos), 



(cos) ~ 



1 



q 2—1 

'H„2 COS) 



(H p cos) 



/? 3 



HJ cos 



(cos) 



(H n 



{ (gp cos) , „, 

+ ^ (cos), — \ h p 2 )* j - 2 ( —^r. — \ p>* 



p{ (H p ), (Hp cos) . _ ^ v ; 



/9 



(cos) ^ 



q,q' i=l j=l 



(if p cos) 4 
(cos) ^ 

i i 



0T o 



\ E E + ? E E E E ^ +q ' r ^ q pq 

q 1=1 q q' 1=1 j=l 

III 



'PqPq'Pq" 



q,q',q" z— 1 j — 1 fc— 1 



(44) 



where in the brackets we have noted for economy of space 
cos -j= Ek (7y£0 k + ?7k^k) = cos • O ne can easu y nn d 



(45) 
i ''k''-k 



where 77k = + i?7 k and the sum runs over the whole 1st 
Brillouin zone except the point k — 0. The quantities of 

the form (<? kl ■ • ■ k „ cos E k + Vlfii)) ^ in 61 

can be obtained using the property: 



\ 2l / di] 



ki 



(46) 



where the notation 



is used. 



Now, taking derivatives step by step and applying the 

o, we obtain: 



obvious equality J^i 



0*0* cos^=E k (« + «) 

cos^E k (« + ^ k ) 
rikr)k> 



(47) 



1 



J k +k',0 



{2c 2 f3J) 2 N 7k7k - 2c 2 /3J 7k 



0k A A s 0k, cos Ek + « 

= (cos^E k (« + «))^ 

1 77 kl 77k 2 ?7k3?7k 4 
(2c 2 /3J) 4 iV 2 7 kl 7k 2 7k 3 7k 4 

[ ^ki+k 2 ,0?7k 3 ?7k 4 + ^ki+k 3 .0'7k 2 '7k 4 

7k 2 7k 3 7k 4 

^ka+ka.O^ki^k, 



(48) 



(2c2/3J) 3 ^V Tk 2 7k 3 7k 4 

^ki+k 4 ,o^7k 2 ?7k3 
7k 2 7k 3 7k 4 

^k 2 +k 4 ,o^7ki?7k3 
7ki7k 2 7k 3 

1 / 5 kl+k2; o5k 3 +k 4 ,0 



(2c 2 /3J) 2 



7kj7k 2 



7ki7k 2 7k 4 

7k!7k 2 7k 3 y 

^ki+k 3 ,0^k 2 +k 4 ,0 

7ki7k 2 

^ki+k 4 ,0^k 2 +k 3 ,0 \ 
7k!7k 2 / 



and 



0kAAAAA 6 cos ^ Ek + «) 
= (cos^E k (« + «)X 

(2c 2 pj) 6 N 3 7 kl 7k 2 7k 3 7k 4 7k 5 7k 6 
1 / <5 kl+ k 2 ,o^k 3 ?7k 4 ?7k 5 »7k 6 



(49) 



(2c 2 (3JfN 2 y 7k 2 7k 3 7k 4 7k s 7k 6 
1 / <5 kl +k 2 ,o<5k 3 +k 4 ,o?7k5 7 7k 6 



(2c 2 (3JYN y 7 kl 7k 3 7k 5 7k 6 

1 / <5ki+k 2 ,0^k 3 +k 4 ,0<5kB+k6,0 



(2c 2 /3J) 3 



7k!7k 3 7k 5 



+ 
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where the sums in the brackets run over all possible ways 

of choosing the pairs of k's in the delta-symbols. Averages 

(0kM„, (fkA a 0k3 0k4>, and (6Mk 2 6»k3 0k 4 0k 5 0k 6 >* can S 2 (R,N) = c 2 

be obtained from (|47]) -(|49 f putting all the 's equal to 

zero. 

Substituting these results in (|44|) and applying (f20|) we 
arrive at the expression: 



k,k — .9k, -k) 



k#0 



2 1-^ 



2 qR 



a 2 N 



q#0 



,cos^E k («+«; 



cos^Ek(«+«: 



(50) S 3 (R,N) 



1 -c 1 
4c 3 2V 2 



y^gk.k'ffk'.k 



k,k' 



7k 



1 /3J 

1 - 3c + 2c 2 
2c 1 



- T 



22V 3 



ffk,k'5k',k"5k",k- 



k,k',k" 

with 5k, k' given by (TT5")) . We have neglected terms that 
vanish in the thermodynamic limit and terms containing 
higher powers of l/(/3J), since we consider low tempera- 
tures. 



C3 



22V 



sin 



2 qR 



a 2jV~ |q| 



a 2 2V ^ |q 

q#0 |M 



&(2V) 



1 



52 + 2~/v E (- 9k > k ~ 5k >- k ) 

k^O 



2 1 



1 



a 2 2V ^ Iql ' 

q#0 



Appendix B 

We are interested in the asymptotic behaviour of the sums: 



q#0 



9d> 



7q 



fl 1 — q,q'flq',q"<?q",q~ 



7q 



q,q',q" 



q#0 



S 2 (2V) = ^^Jq,q'3,'4 , 

q,q' 

= W E .9-q,q'.9q'.q".9q".q^ 



q,q',q" 



with <7q, q ' given by Eq. (|15[) . when R — > oo and 2V — > oo. 

The singularity in J- in the point 9 = defines the 
asymptotic behaviour of the sums. Thus we will have the 
same asymptotic behaviour after expanding the expres- 
sions in the sums for small q's: 



SiCR,JV) = ci 



2 1 v-^ sin 



E 



2 qR 



a 2 N |q| 
q#o IMI 



S 3 (N) 



C3 



2 ^ E E ^ 

2 1^1 
~n2~N ^ T77T ' 



.9k ,-k') 



a 2 2V '—f |q 

q#0 IM ' 

c ii c 2, C3, cx, C2, £3 are constants. 
Numerical calculation gives 



--; E (- 9k > k - 5k,-k) 



0.73 



k^O 



22V ^ ^ 

k#0 k'#0 



.9k,k' (,9k, k' - .9k,-k') ~ 0.27 



To get the asymptotic behaviour of the sums 

civ^ 3R 



'e 1 



111 ' a2iV qt0 lql 

we replace sums over the 1st Brillouin zone in the ther- 
modynamic limit with integrals over continuous variables 
q x , q y , according to the formula 



E 



Na 2 

(2tt) 2 



dq x / dq y + o(N x ) 



We take into account the absence of the terms with 9 = 
in the sums cutting out from the continuous domains of 
integration spaces around the points 9 = with area equal 

tn ^ ■ 
t0 NaT- 
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Fig. 6. The cut out of the point q = 



its argument. In the integral over the rest of the domain, 
e, , we substitute sin 2 with its mean value 1/2, this 



cannot change the asymptotic behaviour for a large R. 
Thus we have simply integrablc functions now and one 
easily finds: 



2 1 

a 2 " AT 



E 



2 qR 



q#0 



1 

2^2 



/■27T 

xdx I dip cos 2 ip 

2£. Jo 



+ 



4tt 2 



dx 



2 71 



dip = c 



fin- 

2tt a 



(53) 



U 



i ' aVN I \aVN' a. j 
q y € l^a'^^Tw ) U ( ^7w' a, 



We neglected here terms containing small values (R/a) 2 /N 
and 1/iV and collected terms with the fixed parameter e 
in the constant C . 
Finally we have: 



We should stress here that the exact value of this area 
is not important for the asymptotic behaviour which is 
the point of our interest, it only must be proportional to 
1/JV. 

After passing to polar coordinates q e ( 2 ^rjj , 
and tp e (0, 27r) one can write 



Si(R,N) 
S 2 (R,N) 
S 3 (R,N) 



= c> I + — In B. 

1 2vr a 

= <L + 0.73— In £ 

2 2tt a 

= c', -0.27— ln^ 

3 2tt a 



Note, that the above estimates hold for finite N as well. 



2 1 



E 



a 2 N ^— ' |q| 

q#0 lHl 



sin 2 ^ i 



2tt 2 



2V¥ 



\/]V 



2T sin 2 ^Rmscp 
d<7 / 



and 



2 1^1 
a 2 TV 2-" |q| 



■ dq 



dip 



(51) 
(52) 



The integral in ||52J) gives ^rhiA. So 



Si{N) = cj + — In A , 

S 2 (A) = c 2 + 0.73^1nA , 

S 3 (A) = c 3 - 0.27-!- In A . 

Now we are interested in the behaviour of the integral 
in (|51[) in the asymptotic case: A — > oo, i? — > oo . After 
change of variables, ^ — ► x, we split the integrals over x 
in two parts: 

„_. R\Ztt 



dx 



dx 



dx 



It is reasonable to assume that the ratio ^H=r is small 

VN 

in the thermodynamical limit, then we choose e small 
enough to change sin(xcos</?) in the limits (~^> e ) by 
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